Link del repositorio de github aquí
#cuando te metes
load("OUTPUT/objetos_actuales.Rdata")
#antes de salir
save.image(file="objetos_actuales.Rdata")
library(tidyjson)
library(rjson)
library(tidyverse)
library(dplyr)
library(ggplot2)
library(lubridate)
library(readr)
library(plotly)
En nuestro seminario vamos a estudiar las altas hospitalarias por ictus, en los diferentes hospitales de Castilla y León durante cuatro años, y como estas, están relacionadas con distintos factores de riesgo.
Lo primero de todo, explicar que los factores de riesgo son características o circunstancias que atentan contra el equilibrio, contra la salud, y que causan enfermedades o muerte. [senado1999factores]
Durante el desarrollo de nuestro seminario vamos a trabajar con datos de dos factores de riesgo, los cuales son: la calidad del aire y valorar si los centros de desintoxicación los podemos considerar como factor de riesgo o como factor de prevención de sufrir un ictus.
La calidad del aire es quizas un factor de riesgo del ictus, quizás poco conocido, pero estudios confirman que “Niveles más elevados de NO el día de ingreso conlleva un aumento del riesgo de mortalidad” [suarez2024influencia]
Objetivos principales:
1.Estudiar la cantidad de altas hospitalarias por ictus que hay en cada provincia de Castilla y León y ver como varían dependiendo del sexo y de la edad.
2.Analizar si a peor calidad del aire hay más riesgo de sufrir un ictus.
3.Estudiar la relación entre centros de ayuda a la desintoxicación por provincia y el número de ictus.
Los datos han sido recabados de Datos abiertos del gobierno de españa que hasta el momento contiene cerca de 90000 conjuntos de datos, distribuidos en un amplio número de temáticas.También supera las 550000 distribuciones, Además de 313 iniciativas de datos abiertos.
Descripción:
Contamos con 1 set de datos en tipo .csv. Este formato, (Comma Separated Values) son archivos de texto que en pricipio deberían estar con los caracteres separados por comas, aunque en algunos de nuestros datos se usa punto y coma en lugar de la coma.
Este es nuestro archivo .csv.:
También contamos con dos sets de datos .json (JavaScript Object Notation) es un formato para almacenar e intercambiar datos.
Estos son nuestros archivos .json.:
-centros_servicios_saisde Este último proviene del Instituto Nacional de Estadística (INE).
altas_hospitalarias_ictus <- fromJSON(file="INPUT/data/altas_hospitalarias_ictus.json")
#altas_hospitalarias_ictus
calidad_aire <- read_csv("INPUT/data/calidad_aire.csv")
#calidad_aire
centros_servicios_saisde <- fromJSON(file = "INPUT/data/centros_servicios_saisde.json")
#centros_servicios_saisde
Importamos la función que cuenta el número de ictus por provincia:
Está función se creo antes de conocer tuberías, la función group_by y summarise con las cuales podemos obetener el mismo resultado pero lo dejaremos a modo comparativo, para ver las ventajas y eficiencia de trabajar con estas estructuras.
source("INPUT/functions/ictus_por_provincia.R")
Aquí vemos el ejemplo trabajando con dichas metodologías
prov_e_ictus <- Provincias%>%
group_by(Provincia)%>%
summarise(numero_ictus=n())
Para resolver el primer objetivo vamos a:
altas_hospitalarias_ictus %>%
spread_all() %>%
gather_object %>%
json_types %>%
count(name, type)
## # A tibble: 4 × 3
## name type n
## <chr> <fct> <int>
## 1 datasetid string 20090
## 2 fields object 20090
## 3 record_timestamp string 20090
## 4 recordid string 20090
Vemos que en fields tenemos un objeto con los datos que queremos extraer, por lo tanto hacemos:
altas_hospitalarias_ictus %>%
enter_object(fields) %>%
gather_object %>%
spread_all %>%
json_types%>%
count(name, type)
## # A tibble: 11 × 3
## name type n
## <chr> <fct> <int>
## 1 ambito_de_procedencia string 20040
## 2 area string 20090
## 3 dia_de_la_semana_en_la_fecha_del_alta string 20090
## 4 dia_de_la_semana_en_la_fecha_del_ingreso string 20090
## 5 edad number 20090
## 6 fecha_de_alta string 20090
## 7 fecha_de_ingreso string 20090
## 8 hospital string 20090
## 9 provincia string 20090
## 10 sexo string 20090
## 11 zona_basica_de_salud_del_paciente string 20038
Dejamos la tabla entera tabulada, para trabajar a partir de ahora con este objeto:
altas_hosp_ictus_tab <- altas_hospitalarias_ictus%>%
spread_all()%>%
enter_object("fields")%>%
spread_all()
#altas_hosp_ictus_tab
Obtener el número de casos por provincias para empezar a relacionarlo con los siguientes objetivos.
Provincias <- altas_hosp_ictus_tab%>%
rename(Provincia=fields.provincia)%>%
select(Provincia)
Provincias$Provincia <- ifelse(Provincias$Provincia == "Ávila", "Avila", Provincias$Provincia)
#Provincias
Obtenemos el número de ictus por provincia
prov_e_ictus <- Provincias%>%
group_by(Provincia)%>%
summarise(numero_ictus=n())
Calculamos el número de personas que han sufrido un ictus para cada edad, para ello obtenemos el atributo edad de cada una de las personas. Además, lo haremos de dos maneras:
edades_total <- altas_hosp_ictus_tab %>%
rename(Edad = edad) %>%
select(Edad)
#edades_total
Importamos función que cuenta el total para cada una de las edades:
Se trata de una función programada desde un punto de vista de programación, más adelante mostraremos como hacerlo desde un enfoque en el manejo de datos como se requiere en la asignatura y por ende, la gran mayoría de código quedará implementado de la segunda manera.
source("INPUT/functions/total_por_edad.R")
resultados_edad <- total_por_edad(edades_total)
resultados_edad
## Edad Conteo
## 1 1 3
## 2 2 0
## 3 3 0
## 4 4 0
## 5 5 1
## 6 6 0
## 7 7 1
## 8 8 3
## 9 9 5
## 10 10 0
## 11 11 0
## 12 12 1
## 13 13 2
## 14 14 4
## 15 15 2
## 16 16 1
## 17 17 2
## 18 18 1
## 19 19 2
## 20 20 3
## 21 21 6
## 22 22 3
## 23 23 2
## 24 24 2
## 25 25 6
## 26 26 2
## 27 27 4
## 28 28 9
## 29 29 5
## 30 30 7
## 31 31 12
## 32 32 8
## 33 33 10
## 34 34 9
## 35 35 15
## 36 36 19
## 37 37 21
## 38 38 28
## 39 39 20
## 40 40 33
## 41 41 36
## 42 42 53
## 43 43 46
## 44 44 53
## 45 45 57
## 46 46 75
## 47 47 79
## 48 48 102
## 49 49 97
## 50 50 106
## 51 51 111
## 52 52 141
## 53 53 143
## 54 54 165
## 55 55 190
## 56 56 180
## 57 57 190
## 58 58 201
## 59 59 269
## 60 60 211
## 61 61 282
## 62 62 283
## 63 63 344
## 64 64 338
## 65 65 322
## 66 66 360
## 67 67 393
## 68 68 347
## 69 69 411
## 70 70 412
## 71 71 421
## 72 72 432
## 73 73 456
## 74 74 518
## 75 75 503
## 76 76 614
## 77 77 619
## 78 78 629
## 79 79 587
## 80 80 661
## 81 81 560
## 82 82 612
## 83 83 665
## 84 84 668
## 85 85 729
## 86 86 732
## 87 87 641
## 88 88 711
## 89 89 590
## 90 90 514
## 91 91 456
## 92 92 411
## 93 93 288
## 94 94 248
## 95 95 181
## 96 96 130
## 97 97 85
## 98 98 67
## 99 99 45
## 100 100 22
## 101 101 16
## 102 102 9
## 103 103 8
## 104 104 6
ictus_por_edad <- edades_total%>%
group_by(Edad)%>%
summarise(num_ictus=n())
ictus_por_edad
## # A tibble: 99 × 2
## Edad num_ictus
## <dbl> <int>
## 1 0 7
## 2 1 3
## 3 5 1
## 4 7 1
## 5 8 3
## 6 9 5
## 7 12 1
## 8 13 2
## 9 14 4
## 10 15 2
## # ℹ 89 more rows
Podemos observar que de esta forma obtenemos menos filas, que de la primera forma, esto se debe a que en esta segunda tabla nos quita las personas con las edades que no han sufrido un ictus.
Regresión polinómica:
regresion_edad_ictus <- lm(num_ictus~poly(Edad,3),data=ictus_por_edad)
summary(regresion_edad_ictus)
##
## Call:
## lm(formula = num_ictus ~ poly(Edad, 3), data = ictus_por_edad)
##
## Residuals:
## Min 1Q Median 3Q Max
## -227.10 -81.20 -21.57 79.80 279.25
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 202.93 10.82 18.759 < 2e-16 ***
## poly(Edad, 3)1 1449.03 107.64 13.462 < 2e-16 ***
## poly(Edad, 3)2 -527.86 107.64 -4.904 3.87e-06 ***
## poly(Edad, 3)3 -1368.52 107.64 -12.714 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 107.6 on 95 degrees of freedom
## Multiple R-squared: 0.7943, Adjusted R-squared: 0.7879
## F-statistic: 122.3 on 3 and 95 DF, p-value: < 2.2e-16
Hemos decidido realizar una regresión polinómica ya que una lineal no se ajusta exactamente lo que queremos y nos daba un R^2 muy bajo, en cambio la polinómica nos da un R^2 de 0.7994 lo cual nos indica que no es un mal modelo. Esto sucede porque en la regresión lineal debe crecer linealmente hasta el infinito pero a ciertas edades es biologicamente imposible que aumente ya que los humanos tenemos una vida limitada. En los resultados del modelor polinomial vemos que la edad es un factor estadísticamente siginificativo.
Gráfico que explica el modelo:
ggplot(data=ictus_por_edad, aes(x=Edad,y=num_ictus))+
geom_point()+
geom_smooth()+
labs(x="Edad (años)",
y="Número de ictus",
title="Regresión polinómica edad-ictus")
Se ve claramente como a edades bajas el número deictus es muy reducido pero a partir de la barrera de los 45-50 empieza a subir drásticamente hasta los 75-88 que empieza a disminuir ya que a estas edades comienza a disminuir la tasa de supervivencia ya que las personas fallecen por otras patologías o por causas naturales.
Ahora vamos a obtener el número de hombres y mujeres para ver si el sexo influye en la cantidad de ictus
Seleccionamos el sexo:
sexototal <- altas_hosp_ictus_tab %>%
rename(Sexo = sexo) %>%
select(Sexo)
total_sexo <- sexototal$Sexo
Llamamos a la función que nos cuenta el número de hombres y de mujeres que han sufrido un ictus:
source("INPUT/functions/total_por_sexo.R")
resultados_sexo <- total_por_sexo(total_sexo)
print(resultados_sexo)
## Sexo Conteo
## 1 Hombre 11124
## 2 Mujer 8966
ictus_por_sexo<- sexototal%>%
group_by(Sexo)%>%
summarise(num_ictus=n())%>%
arrange(desc(Sexo))
ictus_por_sexo
## # A tibble: 2 × 2
## Sexo num_ictus
## <chr> <int>
## 1 Mujer 8966
## 2 Hombre 11124
Un total de 11124 hombres sufrieron un ictus y el número de mujeres fue de 8966.
Gráfico de barras:
Mostramos de manera visual con un gráfico de barras sencillo la cantidad de hombres y mujeres que han sufrido un ictus:
ictus_por_sexo%>%
ggplot(mapping=aes(x=Sexo, y=num_ictus))+
geom_bar(stat="identity", aes(fill=Sexo))+
scale_fill_manual(values=c("Hombre"="blue", "Mujer"="deeppink"))+
labs(x="Sexo",y="Numero de ictus",title="Cantidad de personas (por sexo) que sufrieron un ictus")
ANOVA:
Realizamos un ANOVA para ver si el sexo influye en el riesgo de sufrir un ictus.
Para ello necesitamos conocer la media de la edad y su desviación estandar.
media_desviacion <- altas_hosp_ictus_tab %>%
group_by(sexo) %>%
summarise(
media = mean(edad, na.rm = TRUE),
desviacion_estandar = sd(edad, na.rm = TRUE)
)
media_desviacion
## # A tibble: 2 × 3
## sexo media desviacion_estandar
## <chr> <dbl> <dbl>
## 1 Hombre 73.3 13.0
## 2 Mujer 78.4 13.2
anova_resultado <- aov(edades_total$Edad ~ sexototal$Sexo, data = altas_hosp_ictus_tab)
summary(anova_resultado)
## Df Sum Sq Mean Sq F value Pr(>F)
## sexototal$Sexo 1 131314 131314 765.3 <2e-16 ***
## Residuals 20088 3446613 172
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Análisis estadístico:
Tomamos como Test Hipótesis:
H0: Medias iguales –> p-valor>0.05
H1: Medias distintas –>p-valor<0.05
Tras realizar el anova(analisis de la varianza) con la edad como variable numérica y el sexo como variable categórica vemos que el p-valor es muy bajo rechazamos la H0, lo que nos indica que el sexo es un factor estadisticamente significativo para los ictus.
Calculamos la media de los factores más destacados del aire, agrupados por fecha y provincia.
media_factores_aire <-calidad_aire%>%
rename(provincia=Provincia)%>%
mutate(Fechas = my(Fecha))%>%
mutate(mes= month(Fechas))%>%
mutate(anyo= year(Fechas))%>%
group_by(Fechas, provincia,mes,anyo)%>%
summarise(
NO=mean(`NO (ug/m3)`, na.rm = TRUE),
O3=mean(`O3 (ug/m3)`, na.rm = TRUE),
PM25=mean(`PM25 (ug/m3)`, na.rm = TRUE)
)%>%
arrange(Fechas)
#media_factores_aire
Necesitamos cambiar el formato de la fecha, para así poder trabajar y agrupar datos que contienen fechas, ya que el paquete lubridate permite hacer lo que necesitamos con las fechas, en formato fecha.
altas_hosp_ictus_fecha_bien <- altas_hosp_ictus_tab %>%
mutate(Fecha = ymd_hms(fields.fecha_de_ingreso))%>%
mutate(mes= month(fields.fecha_de_ingreso))%>%
mutate(anyo = year(fields.fecha_de_ingreso))
altas_hosp_ictus_fecha_bien$provincia <- ifelse(altas_hosp_ictus_fecha_bien$provincia == "Ávila", "Avila",altas_hosp_ictus_fecha_bien$provincia)
#altas_hosp_ictus_fecha_bien
Juntamos ambas tablas con sus atributos comunes que nos son de interes (mes, año y provincia)
aire_con_ictus <- full_join(altas_hosp_ictus_fecha_bien, media_factores_aire, by = c("mes","anyo", "provincia"))
aire_con_ictus_resumen <- aire_con_ictus%>%
rename(ambito=fields.ambito_de_procedencia,dia_semana=fields.dia_de_la_semana_en_la_fecha_del_ingreso)%>%
select(Fechas,mes,anyo,edad,sexo,dia_semana,provincia,NO,O3,PM25,ambito)
#aire_con_ictus_resumen
Calculamos la media de cada una de las sustancias (NO, O3 y PM25) agrupadas por provincia, año y mes. Además, contamos el número de ictus que hay para cada uno de ellos.
calcula_media_sustancias_por_provincia_y_meses <-aire_con_ictus_resumen %>%
group_by(provincia, anyo, mes) %>%
summarize(num_ictus=n(),
avg_NO = mean(NO, na.rm = TRUE),
avg_O3=mean(O3, na.rm = TRUE),
avg_PM25=mean(PM25, na.rm = TRUE))
calcula_media_sustancias_por_provincia_y_meses #AÑADIR library(DT)
## # A tibble: 459 × 7
## # Groups: provincia, anyo [45]
## provincia anyo mes num_ictus avg_NO avg_O3 avg_PM25
## <chr> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
## 1 Avila 2019 10 1 2 52 NaN
## 2 Avila 2019 11 1 2 60 NaN
## 3 Avila 2019 12 6 6 43 NaN
## 4 Avila 2020 1 24 5 46 NaN
## 5 Avila 2020 2 17 5 42 NaN
## 6 Avila 2020 3 29 3 66 NaN
## 7 Avila 2020 4 19 2 65 NaN
## 8 Avila 2020 5 20 2 69 NaN
## 9 Avila 2020 6 37 3 72 NaN
## 10 Avila 2020 7 24 2 89 NaN
## # ℹ 449 more rows
Estudios confirman que “Niveles más elevados de NO el día de ingreso conlleva un aumento del riesgo de mortalidad debido a ictus” El gráfico que se muestra a continuación nos sirve para ver sí cuando la concentración de NO es mayor hay un mayor número de ictus.
calcula_media_sustancias_por_provincia_y_meses%>%
ggplot(aes(x=avg_NO, y=num_ictus))+
geom_point(aes(colour = factor(provincia)))+
geom_smooth(method = "lm", aes(colour = factor(provincia)))+
labs(x="Media concentración NO (ug/m3)", y="Número de ictus", title="Gráfico de dispersión del NO", colour="Provincias")
Podemos observar que las lineas para cada provincia según el articulo(ref) deberían salir con pendiente positiva pero no es así. Creemos que este resultado se debe a un pequeño matiz entre el artículo y nuestros datos. Dicho matiz se trata de que en el artículo nos indica que concentraciones altas de NO aumentaban el riesgo de mortalidad mientras que en nuestros datos los pacientes que tenemos han recibido el alta. Así que como no tenemos número de fallecidos no lo podemos contrastar realmente con el artículo.
El ozono es uno de los principales contaminantes del aire. Tal y como afirma el profesor Shaowei Wu, de la universidad Xi´an Jiaotong (China): “estudios previos han sugerido que la contaminación por ozono daña el corazón y los vasos sanguíneos”.
Por ello vamos a realizar un gráfico de dispersión en el que veamos como varía el número de ictus dependiendo de como sea la media de O3 en cada provincia.
Gráfico de dispersión del 03:
calcula_media_sustancias_por_provincia_y_meses%>%
ggplot(aes(x=avg_O3, y=num_ictus))+
geom_point(aes(colour = factor(provincia)))+
geom_smooth(method = "lm", aes(colour = factor(provincia)))+
labs(x="Media concentración O3 (ug/m3)", y="Número de ictus", title="Gráfico de dispersión del O3", colour="Provincias")
En este caso, sí que podemos observar una pendiente positiva a mayor concentración de ozono por lo que podemos concluir que si es un factor que aumenta el numero de ictus.
Las hospitalizaciones aumentan en un día determinado entre un 0,5% y un 1% por cada aumento de 10 mg/m3 de PM2.5. PM2.5: Hace referencia al material particulado respirable presente en la atmósfera. Se puede encontrar más información sobre PM2.5 aqui
Por ello vamos a realizar un gráfico de dispersión para ver si cuando los niveles de PM2.5 son más elevados hay mayor número de ictus.
Gráfico de dispersión del PM25:
calcula_media_sustancias_por_provincia_y_meses%>%
ggplot(aes(x=avg_PM25, y=num_ictus))+
geom_point(aes(colour = factor(provincia)))+
geom_smooth(method = "lm", aes(colour = factor(provincia)))+
labs(x="Media concentración PM2,5 (ug/m3)", y="Número de ictus", title="Gráfico de dispersión del PM2,5", colour="Provincias")
En este caso tenemos provincias que no tienen valores de PM2.5 (las que no tienen linea en el gráfico) pero en las que lo tenemos en todas se cumple que a mayor concentración mayor número de ictus.
Gráfico de barras de todos los gases:
Vamos a graficar todos los gases en una misma imagen, para ello, necesitamos modificar nuesta tabla aire_con_ictus_resumen añadiendo una columna llamada gases que contenga el nombre de los gases que hemos seleccionado y una columna media de concentración con los valores de dicha concentración:
tabla_con_columna_con_todos_gases <- aire_con_ictus_resumen%>%
select(mes,anyo,provincia,NO,O3,PM25)%>%
pivot_longer(names_to="gases",values_to="media_de_concentracion",cols=c(NO,O3,PM25))
tabla_con_columna_con_todos_gases
## # A tibble: 60,306 × 5
## mes anyo provincia gases media_de_concentracion
## <dbl> <dbl> <chr> <chr> <dbl>
## 1 12 2023 Avila NO 4
## 2 12 2023 Avila O3 38
## 3 12 2023 Avila PM25 NaN
## 4 12 2023 Valladolid NO 12.5
## 5 12 2023 Valladolid O3 31
## 6 12 2023 Valladolid PM25 10
## 7 12 2023 León NO 6.88
## 8 12 2023 León O3 34.5
## 9 12 2023 León PM25 9.5
## 10 12 2023 Salamanca NO 6.67
## # ℹ 60,296 more rows
Faceta:
Realizamos una faceta para mostrar lo explicado anteriormente:
tabla_con_columna_con_todos_gases %>%
ggplot(aes(x=anyo, y=media_de_concentracion))+
geom_point(aes(colour = factor(provincia)))+
geom_smooth(method = "lm", aes(colour = factor(provincia)))+
facet_wrap(facets=vars(gases), nrow=1)+
labs(x="Año",y="Media de concentración", title="Faceta de las sustancias del aire", colour="Provincias")
Ahora vamos a realizar 3 gráficos, uno para cada sustancia
Gráfico de barras de la media de NO:
Gráfico que muestra la media de NO (en una escala de colores) por provincia y mes:
calcula_media_sustancias_por_provincia_y_meses %>%
ggplot(aes(x = reorder(factor(provincia),factor(anyo)), y = num_ictus)) +
geom_bar(stat = "identity", aes(fill = avg_NO)) +
scale_fill_gradient(low = "turquoise", high = "blue" , space = "Lab", na.value = "grey50", guide = "colourbar") +
labs(x = "", y = "Número ictus", fill = "Media de concentración de NO") +
theme_classic()
En esta gráfica aparecen representados los números de ictus por
provincia y la concentración de NO, la cual aparece pintada más oscura
cuando es mayor y más clara cuando es menor.
p <- calcula_media_sustancias_por_provincia_y_meses %>%
ggplot(aes(x = reorder(factor(provincia),factor(anyo)), y = num_ictus)) +
geom_bar(stat = "identity",aes(fill = avg_NO)) +
scale_fill_gradient(low = "turquoise", high = "blue", space = "Lab", na.value = "grey50", guide = "colourbar") +
labs(x = "Provincia y Año", y = "Número de ictus", fill = "Media de concentración de NO") +
theme_classic() +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
p
# Convierte el gráfico ggplot a plotly
ggplotly(p, tooltip = c("x", "y", "fill"))
Esta gráfica es una versión mejorada de la anterior, en la cual aparecen representados los ictus que ha habido en cada provincia y en cada año, cada barra tiene 12 partes diferenciadas (una por mes) en las de los años 2020 en adelante; en el año 2019 aparecen menos ya que solo contamos con datos de octubre en adelante. Las barras están coloreadas de distintas tonalidades, corriespondiendose las oscuras a mayores concentraciones de NO y las claras a menores concentraciones de NO. Además pasando el cursor por las barras, gracias a la tecnología del paquete plotly, podemos ver podemos ver con que año y provincia se corresponde esa barra, la media de concentración de NO en ese mes y también el número total de ictus que han ocurrido ese mes.
Gráfico de barras de la media de O3:
Gráfico que muestra la media de 03 (en una escala de colores) por provincia y mes:
calcula_media_sustancias_por_provincia_y_meses %>%
distinct(provincia, .keep_all = TRUE) %>%
ggplot(aes(x = reorder(provincia,anyo), y = num_ictus)) +
geom_bar(stat = "identity", aes(fill = avg_O3)) +
scale_fill_gradient(low = "turquoise", high = "blue", space = "Lab", na.value = "grey50", guide = "colourbar") +
labs(x = "", y = "Número ictus", fill = "Media de concentración de 03") +
theme_classic()
En esta gráfica aparecen representados los números de ictus por
provincia y la concentración de O3, la cual aparece pintada más oscura
cuando es mayor y más clara cuando es menor.
p1 <- calcula_media_sustancias_por_provincia_y_meses %>%
ggplot(aes(x = interaction(provincia, anyo), y = num_ictus, fill = avg_O3)) +
geom_bar(stat = "identity") +
scale_fill_gradient(low = "turquoise", high = "blue", space = "Lab", na.value = "grey50", guide = "colourbar") +
labs(x = "Provincia y Año", y = "Número de ictus", fill = "Media de concentración de O3") +
theme_classic() +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
# Convierte el gráfico ggplot a plotly
ggplotly(p1, tooltip = c("x", "y", "fill"))
Esta gráfica es una versión mejorada de la anterior, en la cual aparecen representados los ictus que ha habido en cada provincia y en cada año, cada barra tiene 12 partes diferenciadas (una por mes) en las de los años 2020 en adelante; en el año 2019 aparecen menos ya que solo contamos con datos de octubre en adelante. Las barras están coloreadas de distintas tonalidades, corriespondiendose las oscuras a mayores concentraciones de O3y las claras a menores concentraciones de O3. Además pasando el cursor por las barras, gracias a la tecnología del paquete plotly, podemos ver podemos ver con que año y provincia se corresponde esa barra, la media de concentración de NO en ese mes y también el número total de ictus que han ocurrido ese mes.
Gráfico de barras de la media de PM2,5:
Gráfico que muestra la media de PM2.5 (en una escala de colores) por provincia y mes:
calcula_media_sustancias_por_provincia_y_meses %>%
distinct(provincia, .keep_all = TRUE) %>%
ggplot(aes(x = reorder(provincia,anyo), y = num_ictus)) +
geom_bar(stat = "identity", aes(fill = avg_PM25)) +
scale_fill_gradient(low = "turquoise", high = "blue", space = "Lab", na.value = "grey50", guide = "colourbar") +
labs(x = "", y = "Número ictus", fill = "Media de concentración de PM2.5") +
theme_classic()
En esta gráfica aparecen representados los números de ictus por
provincia y la concentración de O3, la cual aparece pintada más oscura
cuando es mayor y más clara cuando es menor.
p2 <- calcula_media_sustancias_por_provincia_y_meses %>%
ggplot(aes(x = interaction(provincia, anyo), y = num_ictus, fill = avg_PM25)) +
geom_bar(stat = "identity") +
scale_fill_gradient(low = "turquoise", high = "blue", space = "Lab", na.value = "grey50", guide = "colourbar") +
labs(x = "Provincia y Año", y = "Número de ictus", fill = "Media de concentración de PM2,5") +
theme_classic() +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
# Convierte el gráfico ggplot a plotly
ggplotly(p2, tooltip = c("x", "y", "fill"))
Esta gráfica es una versión mejorada de la anterior, en la cual aparecen representados los ictus que ha habido en cada provincia y en cada año, cada barra tiene 12 partes diferenciadas (una por mes) en las de los años 2020 en adelante; en el año 2019 aparecen menos ya que solo contamos con datos de octubre en adelante. Las barras están coloreadas de distintas tonalidades, corriespondiendose las oscuras a mayores concentraciones de O3y las claras a menores concentraciones de O3. Además pasando el cursor por las barras, gracias a la tecnología del paquete plotly, podemos ver podemos ver con que año y provincia se corresponde esa barra, la media de concentración de NO en ese mes y también el número total de ictus que han ocurrido ese mes. Cabe destacar que en este caso en concreto, hay muchos valores NaN, que son los que aparecen de color gris.
Gráfico de barras de los días de ingreso:
Como dato curioso, queremos ver si el día de la semana está relacionado con los ictus
numero_ictus_dia_semana <-aire_con_ictus_resumen %>%
group_by(dia_semana) %>%
summarize(num_ictus=n())
numero_ictus_dia_semana %>%
ggplot(mapping=aes(x = dia_semana, y = num_ictus)) +
geom_bar(stat = "identity", fill="orchid") +
labs(x = "Días de la semana", y="Número de ictus") +
theme_classic()
Por lo leído y descrito en uno de los artículos, niveles altos de NO el día de ingreso conlleva a mas riesgo de mortalidad, al ser nuestros datos de altas hospitalarias por ictus y no de muertes, lo que vamos a intentar ver es si los días que hay un mayor número de ingresos coinciden con niveles altos de NO.
Importamos los datos de los centros de desintoxicación.
Visionamos los datos y vemos que con hacer un spread_all podemos obtener lo que necesitamos para nuestro trabajo.
#spread_all(centros_servicios_saisde)
Centros <- centros_servicios_saisde %>%
spread_all() %>%
mutate(fields.provincia = ifelse(fields.provincia == "Ávila", "Avila", fields.provincia)) %>%
count(provincia = fields.provincia) %>%
as.data.frame()
Centros #UDAR LIBRARY(DT)
## provincia n
## 1 Avila 6
## 2 Burgos 12
## 3 León 18
## 4 Palencia 13
## 5 Salamanca 16
## 6 Segovia 5
## 7 Soria 4
## 8 Valladolid 17
## 9 Zamora 8
ordenadas <- arrange(Centros,Centros$provincia)
ordenadas #UDAR LIBRARY(DT)
## provincia n
## 1 Avila 6
## 2 Burgos 12
## 3 León 18
## 4 Palencia 13
## 5 Salamanca 16
## 6 Segovia 5
## 7 Soria 4
## 8 Valladolid 17
## 9 Zamora 8
Importamos función que imprime los habitantes de castilla y león por provincias:
source("INPUT/functions/poblacion_cyl.R")
poblacion_cyl #UDAR LIBRARY(DT)
## provincia Poblacion_total Poblacion_extranjera
## 1 Avila 159764 12868
## 2 Burgos 357370 33150
## 3 León 448573 23398
## 4 Palencia 157787 8581
## 5 Salamanca 327089 18606
## 6 Segovia 155332 20015
## 7 Soria 89528 10363
## 8 Valladolid 521333 33352
## 9 Zamora 166927 7601
## 10 Castilla y León 2383703 167934
## Porc_pob_extranj_sobre_pob_total Num_munic_pob_extranj_mas_3_por_ciento
## 1 8.05 106
## 2 9.27 220
## 3 5.21 102
## 4 5.43 88
## 5 5.68 130
## 6 12.88 167
## 7 11.57 83
## 8 6.39 150
## 9 4.55 110
## 10 7.04 1156
#str(poblacion_cyl)
Para tener una referencia entre los centros de desintoxicación por provincia con el número de habitantes de esa provincia calcularemos cuál es el número de habitantes por centro; con el obetivo de ver realmente si hay más ictus en aquellas provincias donde el número de habitantes por centro es menor y el consumo de drogas es un factor significativo en los ictus o por el contrario no lo es.
pob_centros <- full_join(poblacion_cyl,ordenadas)
tabla <- pob_centros%>%
select(provincia,Poblacion_total,n)
source("INPUT/functions/divide.R")
hab_centro <- divide(tabla$Poblacion_total,tabla$n)
#faltaria indicar que con que provincia se corresponden esos datos.
tabla$hab_centro <- hab_centro
resultado <- tabla %>%
select(provincia, hab_centro)
resultado #UDAR LIBRARY(DT)
## provincia hab_centro
## 1 Avila 26627.33
## 2 Burgos 29780.83
## 3 León 24920.72
## 4 Palencia 12137.46
## 5 Salamanca 20443.06
## 6 Segovia 31066.40
## 7 Soria 22382.00
## 8 Valladolid 30666.65
## 9 Zamora 20865.88
## 10 Castilla y León NA
Lo lógico es esperar que en las provincias que cuentan con un menor número de habitantes, haya menor número de ictus y menor número de centros de desintoxicación.
numero_ictus <- ictus_por_provincia(Provincias)
df <- as.data.frame(numero_ictus)
pob_prov_y_total <- poblacion_cyl%>%
select(provincia, Poblacion_total)
pob_prov_y_total
## provincia Poblacion_total
## 1 Avila 159764
## 2 Burgos 357370
## 3 León 448573
## 4 Palencia 157787
## 5 Salamanca 327089
## 6 Segovia 155332
## 7 Soria 89528
## 8 Valladolid 521333
## 9 Zamora 166927
## 10 Castilla y León 2383703
pob_cyl_sin_ult_fila <- pob_prov_y_total [1:9,]
pob_cyl_sin_ult_fila
## provincia Poblacion_total
## 1 Avila 159764
## 2 Burgos 357370
## 3 León 448573
## 4 Palencia 157787
## 5 Salamanca 327089
## 6 Segovia 155332
## 7 Soria 89528
## 8 Valladolid 521333
## 9 Zamora 166927
arrange(.data=ordenadas,(ordenadas$provincia))
## provincia n
## 1 Avila 6
## 2 Burgos 12
## 3 León 18
## 4 Palencia 13
## 5 Salamanca 16
## 6 Segovia 5
## 7 Soria 4
## 8 Valladolid 17
## 9 Zamora 8
datos <- data.frame(pob=pob_cyl_sin_ult_fila,
centros= arrange(.data=ordenadas,(ordenadas$provincia)),
num_ictus=df$numero_ictus)
str(datos)
## 'data.frame': 9 obs. of 5 variables:
## $ pob.provincia : chr "Avila" "Burgos" "León" "Palencia" ...
## $ pob.Poblacion_total: num 159764 357370 448573 157787 327089 ...
## $ centros.provincia : chr "Avila" "Burgos" "León" "Palencia" ...
## $ centros.n : int 6 12 18 13 16 5 4 17 8
## $ num_ictus : num 1178 2788 3513 1064 3743 ...
datos
## pob.provincia pob.Poblacion_total centros.provincia centros.n num_ictus
## 1 Avila 159764 Avila 6 1178
## 2 Burgos 357370 Burgos 12 2788
## 3 León 448573 León 18 3513
## 4 Palencia 157787 Palencia 13 1064
## 5 Salamanca 327089 Salamanca 16 3743
## 6 Segovia 155332 Segovia 5 1142
## 7 Soria 89528 Soria 4 601
## 8 Valladolid 521333 Valladolid 17 4575
## 9 Zamora 166927 Zamora 8 1486
Gráfico de barras en el que aparece representado el número de centros por provincia:
datos%>%
ggplot(mapping = aes(x = reorder(pob.provincia, num_ictus), y = centros.n)) +
geom_bar(stat = "identity",fill="blue") +
labs(x = "Provincias", y = "Número de centros") +
theme_classic()
Gráfica de barras en la que se representa el número de ictus en cada provincia:
Cabe destacar que se coloreande un color más fuerte las barras con más centros de desintoxicación.
datos %>%
ggplot(mapping = aes(x = reorder(pob.provincia, num_ictus), y = num_ictus)) +
geom_bar(stat = "identity", aes(fill = centros.n)) +
scale_fill_gradient(low = "peachpuff", high = "red") +
labs(x = "Provincias", y = "Número de ictus", fill = "Centros detox") +
theme_classic()
Gráfico de barras en el que representamos el número de habitantes por centro para cada provincia:
resultado
## provincia hab_centro
## 1 Avila 26627.33
## 2 Burgos 29780.83
## 3 León 24920.72
## 4 Palencia 12137.46
## 5 Salamanca 20443.06
## 6 Segovia 31066.40
## 7 Soria 22382.00
## 8 Valladolid 30666.65
## 9 Zamora 20865.88
## 10 Castilla y León NA
resultado_sin_ult_fila=resultado[1:9,]
str(resultado_sin_ult_fila)
## 'data.frame': 9 obs. of 2 variables:
## $ provincia : chr "Avila" "Burgos" "León" "Palencia" ...
## $ hab_centro: num 26627 29781 24921 12137 20443 ...
ggplot(mapping = aes(x = reorder(datos$pob.provincia, datos$num_ictus), y = resultado_sin_ult_fila$hab_centro)) +
geom_bar(stat = "identity", fill = "green") +
labs(x = "Provincias", y = "Número habitantes por centro") +
theme_classic()
Obtenemos una tabla que nos permite visualizar de manera numérica toda la información estudiada anteriormente:
tabla_centros <- tibble(Provincia=datos$pob.provincia,
Numero_de_ictus=datos$num_ictus, Poblacion_total_cyl=pob_cyl_sin_ult_fila$Poblacion_total,
Numero_centros=datos$centros.n,
Habitantes_por_centro=resultado_sin_ult_fila$hab_centro,
Porcentaje_ictus_por_hab=(Numero_de_ictus/Poblacion_total_cyl)*100)
tabla_centros
## # A tibble: 9 × 6
## Provincia Numero_de_ictus Poblacion_total_cyl Numero_centros
## <chr> <dbl> <dbl> <int>
## 1 Avila 1178 159764 6
## 2 Burgos 2788 357370 12
## 3 León 3513 448573 18
## 4 Palencia 1064 157787 13
## 5 Salamanca 3743 327089 16
## 6 Segovia 1142 155332 5
## 7 Soria 601 89528 4
## 8 Valladolid 4575 521333 17
## 9 Zamora 1486 166927 8
## # ℹ 2 more variables: Habitantes_por_centro <dbl>,
## # Porcentaje_ictus_por_hab <dbl>
# Calcular métricas adicionales
tabla_metricas <- tabla_centros %>%
mutate(tasa_incidencia = (Numero_de_ictus / Poblacion_total_cyl) * 1000,
densidad_centros = (Numero_centros / Poblacion_total_cyl) * 10000)
ggplot(tabla_metricas, aes(x = Numero_centros, y = Numero_de_ictus)) +
geom_point() +
geom_smooth() +
labs(title = "Relación entre numero de centros y numero de ictus",
x = "Centros s",
y = "Número de ictus ")
modelo_regresion_polinomico <- lm(
Numero_de_ictus ~ poly(Numero_centros,4),
data = tabla_metricas
)
summary(modelo_regresion_polinomico)
##
## Call:
## lm(formula = Numero_de_ictus ~ poly(Numero_centros, 4), data = tabla_metricas)
##
## Residuals:
## 1 2 3 4 5 6 7 8
## -229.86 977.06 -333.56 -1180.64 -10.78 -11.33 89.28 586.70
## 9
## 113.14
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2232.2 282.7 7.895 0.00139 **
## poly(Numero_centros, 4)1 3547.4 848.2 4.182 0.01389 *
## poly(Numero_centros, 4)2 633.9 848.2 0.747 0.49637
## poly(Numero_centros, 4)3 183.0 848.2 0.216 0.83972
## poly(Numero_centros, 4)4 -839.9 848.2 -0.990 0.37815
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 848.2 on 4 degrees of freedom
## Multiple R-squared: 0.8267, Adjusted R-squared: 0.6533
## F-statistic: 4.769 on 4 and 4 DF, p-value: 0.07972
https://fhernanb.github.io/libro_regresion/predict.html
tabla_predicciones <- data.frame(tabla_metricas, predicciones = predict(modelo_regresion_polinomico))
# Visualización del impacto
ggplot(tabla_predicciones, aes(x = Habitantes_por_centro, y = tasa_incidencia)) +
geom_point() +
geom_line(aes(y = predicciones), color = "blue") +
labs(title = "Impacto estimado de centros detox en la tasa de ictus",
x = "Habitantes por centro",
y = "Tasa de Ictus por 1,000 habitantes")
Tras la realización de las tres gráficas anteriores para ver como están relacionados el número de ictus de cada provincia con el número de centros de desintoxicación que hay en las mismas, y el número de habitantes por centro de cada una de ellas, hemos llegado a las siguientes conclusiones:
Las poblaciones que cuentan con un menor número de habitantes (Ávila, Palencia, Soria y Zamora) son a su vez las que cuentan con un menor porcentaje de ictus; esto seguramente pueda deberse a que por regla general en las ciudades pequeñas se suelen llevar mejores habitos de vida.
Podemos destacar tres casos en especial:
En Palencia, cuentan con un número elevado de centros de desintoxicación para la población que tiene; lo cual supone que haya pocos habitantes por centro, pudiendo ver así como el porcentaje de ictus es el segundo más pequeño.Por lo tanto podría parecer que la presencia de más centros ayuda a reducir el número de ictus.
Segovia tiene un gran número de habitantes por centro, el más alto de todos, lo cual podría indicar que las personas que acuden van a recibir un “peor” tratamiento ya que el tratamiento se va a tener que repartir entre más pacietes y al ser peor el riesgo de ictus va a ser mayor, sin embargo en esta población no se cumple la hipotésis planteada en la conclusión de Palencia.
Salamanca es la provincia que cuenta con un mayor porcentaje de ictus por habitante, lo cual es destacable ya que no es la que más habitantes tiene y además es la segunda que menos habitantes por centro tiene.
Por lo tanto,visualmente y falta de un estudio estadístico preciso no podemos determinar si el número de centros es un factor estadísticamente significativo para el riesgo de padecer un ictus.
Queda demostrado que tanto la edad como el sexo influyen en el riesgo de sufrir un ictus, esto lo hemos comprobado realizando una regresión polinómica (edad) y un ANOVA (sexo), con sus correspondientes gráficas para hacerlo más visual.
En cuanto a las sustancias del aire, podemos decir que la concentración de NO y la de PM2,5 si que se relacionan con la probabilidad de sufrir un ictus; mientras que la concentración de O3 no podemos llegar a saber si verdaderamente esta relacionada ya que lo que está comprobado es que el hecho de que haya elevadas concentraciones de NO aumenta la mortalidad por ictus en ese día, y nuestros datos no son de muertes sino de altas hospitalarias.
En relación con la cantidad de centros, tampoco podemos determinar si el que haya más numero de centros reduce la probabilidad de sufrir un ictus.
.